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Summary 

- An elliptic grid generation method is presented which generates excellent boundary conforming 
grids in domains in 2D physical space. The method is based on the composition of an algebraic and 
elliptic transformation. The composite mapping obeys the familiar Poisson grid generation system 
with control functions specified by the algebraic transformation. New expressions are given for the 
control functions. Grid orthogonality at the boundary is achieved by modification of the algebraic 
transformation. 

- It is shown that grid generation on a minimal surface in 3D physical space is in fact equivalent 
to grid generation in a domain in 2D physical space. 

- A second elliptic grid generation method is presented which generates excellent boundary con- 
forming grids on smooth surfaces. It is assumed that the surfaces are parametrized and that the 
parametrization is a smooth mapping from a unit square onto the surface. A generated surface 
grid only depends on the shape of the surface and is independent of the parametrization. 

- Concerning surface modeling, it is shown that bicubic Hermite interpolation is an excellent method 
to generate a smooth surface which is passing through a given discrete set of control points. In 
contrast to bicubic spline interpolation, there is extra freedom to model the tangent and twist 
vectors such that spurious oscillations are prevented. 

1 Introduction 

A flow simulation system for the computation of flows about complete aircraft configurations in- 
cluding propulsion aircraft components has been developed at NLR. The system is known as the 
ENFLOW ( Euler/Navier-Stokes FLOW) system [1]. Fig. 4 shows the layout of the system and 
summarizes its use for CFD work. 

Surface modeling of the original aerodynamic input configuration surfaces is done with the commer- 
cial ICEM-CFD software. The subdivision of a three dimensional flow domain into blocks is done 
with the graphical interactive domain modeler ENDOMO. The computation of structured grids in 
the interior of the blocks is done with the graphical interactive grid generator ENG RID. Given a 
multi-block grid, the flow solver ENSOLV computes the solution of the Euler and/or Navier-Stokes 
equations with respect to specified boundary conditions. 

In this paper we focus on surface grid generation. Output of ICEM-CFD is a set of discrete 

1 This investigation was partly performed under contract 01I05N with the Netherlands Agency for Aerospace 
Programs (NIVR). 


mxxam vtm bum tot nura 617 



surfaces. A discrete surface is a two-dimensional array of so-called control points. The complete 
set of discrete surfaces approximates the original aerodynamic input configuration surfaces. 

A discrete surface must be interpolated during grid generation. For this purpose, bicubic Hermite 
interpolation is used to define a smooth surface which is passing through the set of control points 
without introducing spurious oscillations. A C 1 parametrization is constructed which maps a unit 
square onto the interpolated surface. Bicubic Hermite interpolation is discussed in Section 5. 

Multi-block grid generation proceeds from the inside out, starting with the generation of grids in 
edges, followed by the grid generation in faces, and ending with the grid generation in blocks. For 
this reason, a surface grid generation method is needed to generate interior grids in surfaces with 
a prescribed boundary grid point distribution as Dirichlet boundary condition. 

In Section 4 it is shown how elliptic surface grid generation can be used to generate an interior 
surface grid on a parametrized surface with a prescribed boundary grid point distribution. A 
generated surface grid is independent of the parametrization. Thus the interior surface grid will 
only depend on the shape of the surface and the prescribed boundary grid points. 

For surfaces in the interior of a flow domain, often only the boundary shape is defined. For a 
surface with only a prescribed boundary shape and a prescribed boundary grid point distribution, 
it is possible to generate an interior surface grid on a minimal surface. The shape of the surface 
becomes a soap film bounded by the prescribed boundary of the surface. It appears that surface 
grid generation on minimal surfaces is in fact a straightforward extension of grid generation in 2D. 

Grid generation in a domain in 2D is treated in Section 2, the extension to minimal surface grid 
generation is discussed in Section 3. 

Concerning grid generation, the main emphasis lies on the derivation of the elliptic grid generation 
systems. The discretization and solution of the systems of partial differential equations is not 
discussed. We refer to [2] for details about the applied solution methods. Grid generation in 3D 
domains is also discussed in [2]. 


2 2D grid generation 

Consider a simply connected bounded domain V in two dimensional space with Cartesian coor- 
dinates x = ( x,y ) T . Suppose that V is bounded by four edges E\,E 2 ,E 3 ,E 4 . Let (E\,E 2 ) and 
( E 3 , E 4 ) be the two pairs of opposite edges as shown in Fig. 1 . 

Introduce the parameter space V as the unit square in a two dimensional space with Cartesian 
coordinates s = ( s,t) T . Require that the parameters s and t obey: 


• s ee 0 at edge E\ and s = 1 at edge E 2 , ( 1) 

• s is the normalized arc length along edges E 3 and E 4 , (2) 

• t = 0 at edge E 3 and t = 1 at edge E 4 , (3) 

• t is the normalized arc length along edges E\ and E 2 . (4) 
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Computational space C 


Parameter space P 


Domain D 


Figure 1: Transformation from computational (£, tj) space to a domain V in Cartesian ( x,y ) space. 


The mapping s : d D i — * d'P is defined by these requirements. In the interior of T) we require that s 
and t are harmonic functions of x and y, thus obey the Laplace equations: 


As — 

d 2 s 

d 2 s 



dx 2 + dy 2 

— $xx Syy — 0? 

(5) 

At = 

dH 

d 2 t 



dx 2 dy 2 

— ^xx ~l~ ^ yy — 0* 

(6) 


The two Laplace equations A s — 0 and A t = 0, together with the above specified boundary 
conditions, define the mapping s : T> ^ V. Note that this mapping only depends on the shape 
of domain V. By interchanging the dependent and independent variables, a non-linear elliptic 
partial differential equation can be derived for x : V (->• V. This mapping is called the elliptic 
transformation. It is well known that this mapping is differentiable and one-to-one [3]. 

Define the computational space C as the unit square in a two dimensional space with Cartesian 
coordinates £ = (£, r]) T . Assume that a mapping x : dC i-» dT> is prescribed which maps the 
boundary of C one-to-one on the boundary of V. This mapping defines the boundary grid point 
distribution. Assume that 

• £ = 0 at edge E\ and £ = 1 at edge E 2 , 

• 7/ = 0 at edge E 3 and rj = 1 at edge E 4 . 

We wish to construct a mapping i:Ch V which obeys the boundary conditions and which is a 
differentiable one-to-one mapping. Furthermore, we require that the interior grid point distribution 
describes a smooth transition between the prescribed grid point distribution in the four edges. 

A natural mapping x : C >— *■ D exists which obeys these requirements. This mapping will be the 
composition of an algebraic transformation and the elliptic transformation based on the Laplace 
equations. The algebraic transformation is a differentiable one-to-one mapping from computational 
space C onto the parameter space V. The composition of these two mappings defines the interior 
grid point distribution and is a differentiable one-to-one mapping from computational domain C 
onto domain T>. 
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The algebraic transformation is defined as follows. Because x : dC »— > dV is prescribed and x : 
dV dV is defined as described above, it follows that s : dC dV is also defined. 

From the preceding requirements it follows that 

5(0,7;) = 0 , 5(1,7/) = 1 , 5(£,0) = 5£ 3 (f) , s(f, 1) = 5£ 4 (0, (7) 

where the functions se 3 ,se 4 are monotonically increasing, and 

i(£, 0) = 0 , t(t, 1) = 1 , 1(0, 7/) = t El (v) , = tE 2 (V ), (8) 

where the functions are also monotonically increasing. Thus the four functions !#,( 7 /), 

fg 2 ( 7/), 5 £ 3 (£), s e<{£) are defined by the boundary grid point distribution. 

The mapping s : C ^ V is now defined by the following two algebraic equations: 

5 = 5£ 3 (f)(l ~ t) + 5£ 4 (^)1, (9) 

t = *£,(»!)( 1 - 5) + <e 2 (7/)5. (10) 


Eq.(9) implies that a coordinate line £ = constant is mapped to the parameter space V as a straight 
line: s is a linear function of 1, and Eq.(10) implies that a grid line 77 = constant is also mapped 
to as a straight line: t is a linear function of s. For given values of £ and 1 7 , the corresponding 
s and t values are found as the intersection point of the two straight lines. For this reason, the 
system defined by Eqs.(9),(10) is called the “algebraic straight line transformation”. It can be easily 
verified that this system defines a differentiable one-to-one mapping because of the positiveness of 
the Jacobian: s^t v — s v t^ > 0. 

In the remainder of this section, we will derive the set of non-linear elliptic partial differential 
equations which the composite mapping x = x (.?(£)) has to fulfill. First introduce the two covariant 
base vectors 

OX _ _ dx 

fll = at ’ 02 = ait = Xr> ' {11) 

and define the covariant metric tensor components as the inner product of the covariant base vectors 


— (ctj , dj ) , i — ( 1 , 2 ) , j ( 1 , 2 ). 


( 12 ) 


Then the contravariant base vectors a 1 and a 2 are defined according to the rules 


(S \S j ) = 6' j , * = { 1 , 2} , j = { 1, 2), 


(13) 


with Sj the Kronecker symbol. Define the contravariant metric tensor components 


a tJ = (o\a J ) , i = ( 1 , 2 ), j = ( 1 , 2 ), 


(14) 


so that 


j «11 

Y «12 a 22 


a" 

a ' 2 

a ' 2 

a 22 
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and 


a 1 = a n ai + o 12 aj , a 2 = a 12 a 1 + a 22 02. 


(16) 
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Introduce the determinant J 2 of the covariant metric tensor: J 2 = 2 - af 2 . 

Now consider an arbitrary function (f> = 4>{^rf). Then 4> is also defined in domain V and the 
Laplacian of <t> is expressed as 

A <}> = 4> xx + rf* = ±{ (ja u <f>^ + Ja X2 <j> v ^ + (ja} 2 ^ + Ja 22 <j> v ^ | , (17) 

which may be found in every textbook on Tensor Analysis and Differential Geometry (for example 
see [4], page 227). Take as special cases respectively <f> = £ and <f> = T). Then Eq.(17) yields 

A{= + ■ A,) = 7 {(■ / “'*)f + (■'“”)„} ■ (18) 

Thus the Laplacian of 0 can also be expressed as 

A <f> = + 2 a X2 <f> iv + a 22 4> m + Af <^ + ^ <f> n . (19) 

Substitution of respectively (j> = s and <f> = t in this equation yields 

As = a n s^ + 2a X2 s^ v + a 22 s w + A£ + Ai] s v , (20) 

At = a u t^ + 2a 12 t ir) + a 22 t nri + A^t ( + Ajit v . (21) 


Using these equations and the requirement that s and t are harmonic in domain T>, thus As = 0 
and A t = 0, we find the following expressions for the Laplacian of £ and rj 



i xx P n + 2 a X2 P u + a 22 P, 


22 ? 


(22) 


where 



and the matrix T is defined as 




s v 

t v 


(23) 

(24) 


The six coefficients of the vectors P u = ( P } x , P 2 X ) T , P\ 2 = (f^, P 2 2 ) T and P 2 2 = (P^ P 22 V are 
so called control functions. These six control functions are completely defined and easily computed 
for a given algebraic transformation s = s(£). Different and less useful expressions of these control 
functions can also be found in [5, 6]. 


Finally, substitution of cp = x in Eq.(19) yields 


Ax = a u x^ + 2 a x2 x^ n + a 22 x vv + A£ x ^ + At) x n . 


(25) 


Substitution of Eq.(22) into this equation and using the fact that Ax = 0 we arrive at the following 
Poisson grid generation system 


a xx x^ + 2a x2 x{ v + a 22 x vv + 

+ 


(a 11 P/ t + 2 a 12 P 1 1 2 + u 22 P22) 

(a n P 2 + 2o 12 P 1 2 2 + a 22 P 2 2 2 ) *„ = 0. 


(26) 
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Using Eqs.( 12 ),( 15 ), we find the following well known expressions for the contravariant metric 
tensor components: 

J 2 a u = a 2 2 = (%% , J 2 a } 2 = -a 12 = ~(x ( ,x v ) , J 2 a 22 = a u = (%:%. ( 27 ) 

Thus the Poisson grid generation system defined by Eq.( 26 ) can be simplified by multiplication 
with J 2 . Then we obtain: 

022% - 2ai2% + an% + (a n Ph ~ 2ai2%> + <*11^22) H 

+ (&22P11 ~ ^ a \2P\2 + a ii C22) = 0. ( 28 ) 

This equation, together with the expressions for the control functions Pjj given by Eq.( 23 ), forms 
our 2D grid generation system. Grids are computed by solving this quasi-linear system of elliptic 
partial differential equations. 

Orthogonality at boundaries 

Grid orthogonality at boundaries can be achieved as follows. 

Redefine the elliptic transformation x : V ■— ► T> by imposing the following new set of boundary 
conditions for the harmonic functions s and t: 

• s = 0 at edge E\ and s = 1 at edge E2, 

• = 0 along edges £3 and £4, where n is the outward normal direction, 

• t = 0 at edge £3 and t = 1 at edge £4, 

% — 0 along edges E\ and £ 2 , where n is the outward normal direction. 

These new boundary conditions define a new mapping x : V > V. 

The Neumann boundary condition ^ = 0 along edges £3 and £4 imply that a parameter line s = 
constant is a curve in domain V which is orthogonal at those edges. Similarly, a parameter line t = 
constant is a curve in V which is orthogonal at edge E\ and edge £ 2 . 

The algebraic transformation s : C 1— ► V is redefined according to the following two algebraic 
equations: 


S = SE 3 (Z)Ho(t) + ( 29 ) 

< = iEi(v)Po( s ) + t E 2 (v)Pi( s ), ( 30 ) 

where £0 and Hi are cubic Hermite interpolation functions defined in Eq.( 52 ) below. Note that 
£ o ( 0 ) = 1 , ^( 0 ) = 0, Hq( 1 ) = 0, //'( 1) - 0 and Hi( 0 ) = 0 , £{( 0 ) = 0 , H x ( 1 ) = 1 , H[( 1) = 0. 
It follows from Eq.( 29 ) that a coordinate line f = constant is mapped to parameter space £ as a 
cubic curve which is orthogonal at both edge £3 and edge £4 in V . Such a curve in parameter 
space V will thus be mapped by the new elliptic transformation x : V »-► V as a curve which is 
orthogonal at both edge £3 and edge £4 in V . Similar observations can be made for coordinate 

lines 7 ] = constant. Thus the grid will be orthogonal at all four edges in domain V. 
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The composite mapping x : C <-> V still obeys the Poisson grid generation system defined by 
Eq.(28). Thus the same system of elliptic equations can be solved to generate an orthogonal grid 
at the boundary. The only difference is that now s : C <-* V is defined by Eqs.(29),(30) instead of 
Eqs.(9),(10). 

Figs. 5, 6, 7 are demonstrations of the robustness of this elliptic grid generation method. The bound- 
ary grid point distribution is prescribed and the interior grids are obtained by solving Eq.(28). The 
interior grid point distributions were verified to be foldfree by zooming into regions where the grid 
is very dense. 

3 Surface Grid Generation on Minimal Surfaces 

Grid generation on a minimal surface in 3D physical space is in fact equivalent to grid generation 
in a domain in 2D physical space. 

As in the two dimensional case, again consider four curved edges E],E 2 , E 3 , E 4 but now situated 
in the three dimensional physical space with Cartesian coordinates x = (x,y, z) T . Let ( E\,E 2 ) and 
( E 3 , E 4 ) be the two pairs of opposite edges as shown in Fig.2. 

Introduce the parameter space V as the unit square in a two dimensional space with Cartesian 
coordinates s = (s, t) T . Again require that the parameters s and t obey the boundary equations 
specified in Eqs.(l),(2),(3),(4). Furthermore, require that 


As 

= 0, 

(31) 

A t 

- 0 , 

(32) 

H 

= 0, 

(33) 


where A is the Laplace- Beltrami operator for surfaces and H is the mean curvature. 

These three requirements, together with the specified boundary conditions define a unique mapping 
x : V *-* 7 l 3 . The shape of the surface defined by this mapping is a minimal surface (soap film) 
because of the requirement that the mean curvature H is zero. The parametrization of the surface 
is defined by Eqs.(31),(32). Define the minimal surface S = {x(s,£) | (s,£) £ V}. 

Consider a prescribed boundary grid point distribution at the four edges E\, E 2 , £3, E 4 of the 
minimal surface <S. The boundary grid point distribution can be defined as a mapping x : dC dS 
where C is the computational space defined as the unit square in a two dimensional space with 
Cartesian coordinates ( = (£, t]) T . Because x : dC dS is prescribed and x : dV >-> dS is defined 
as described above, it follows that s : dC d~P is also defined. 

In exactly the same way as for the two dimensional case, the mapping s : C t— » 'P is defined by the 
algebraic straight line transformation defined by Eqs.(9),(10). The mapping x : V t-+ S is defined 
by Eqs.(31),(32),(33). The composite mapping x : C S is defined as x - x (s(£)) and describes 
the interior grid point distribution on the minimal surface S. Note that this composite mapping 
will be differentiable and one-to-one. 

We will now show that the set of non-linear elliptic partial differential equations which the composite 
mapping has to fulfill is the same Poisson system as defined by Eq.(28) but with x = ( x,y,z) T 
instead of x = (x,y) T . Thus grid generation on a minimal surface in 3D physical space is in fact 
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Computational space C 


Parameter space P 


Minimal surface S 


Figure 2: Transformation from computational (£,77) space to a minimal surface S in Cartesian 
(x, y, z) space. 


equivalent to grid generation in a domain in 2D physical space. The result that a Poisson system 
of the form as defined by Eq.(28) can be used to compute a grid on a minimal surface can also be 
found as a special application of the formulas derived in [7]. 

For this purpose, introduce the two covariant base vectors 

a x = Xf. , a 2 = x v . (34) 


The two covariant base vectors span the tangent plane of S at the corresponding point P. Define 
the unit surface normal as 

«i A a 2 

n — 7—; — n-, (35) 

II Gi A a 2 || 

where A is the vector product operator. The contravariant base vectors a 1 and a 2 are defined 
according to the rules 

and 


(a l ,a J ) = 6' ] , i = {1,2} , j = {1,2}, 
(a 1 , n) = 0 , ( a 2 , n) = 0. 


(36) 

(37) 


Thus the two contravariant base vectors are also lying in the tangent plane of $ at the corresponding 
point P. Define the covariant metric tensor components by Eq.(12) and the contravariant metric 
tensor components by Eq.(14). Then Eqs.(15),(16) are still valid. Again introduce the determinant 
J 2 of the covariant metric tensor: J 2 = aua 2 2 ~ a ?2* 


Now consider an arbitrary function <f> — <£(£,77). Then <f> is also defined on surface S and the 
Laplace- Beltrami operator of <f> is expressed as 


A <j> =I{ ^Ja 11 ^ + Ja 12 <t> v ^ + (^Ja 12 4>£ + Ja 22 ^^| (38) 

(see [4], page 227). As in the two-dimensional case, substitution of <f> = £ and <fi = 7] into this 
equation yields Eq.(18). Thus the Laplace- Beltrami operator of (j) can also be expressed as defined 
by Eq.(19). Substitution of respectively (j) = s and <f> = t in Eq.(19) and using the requirements ex- 
pressed by Eqs. (31), (32) yields exactly the same expressions for and A77 given by Eqs.(22),(23). 
Finally, substitution of <f> = x in Eq.(19) yields Eq.(25). 
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The Laplace- Beltrami operator applied on x obeys a famous relation expressed by 

A£=2 Hn, (39) 

where the mean curvature H is defined as 

H = ^ (a n % + 2 a n x£ n + a 22 x vn , n) . (40) 

(for example see [8], Theorem 1, page 71). Using the requirement H = 0 yields 

Az = 0. (41) 

Thus Eq.(22) and Eq.(25) with Az = 0 are also valid for minimal surfaces. Following the same 
derivation as given in Section 2, we arrive at exactly the same non-linear system of elliptic partial 
differential equations as expressed by Eq.(28). Thus an interior grid point distribution on a minimal 
surface is found by solving Eq.(28). The only difference compared to the two dimensional case is 
that now z = (z, y , z) T instead of x — (z, y) T . 

Grid orthogonality at boundaries can be obtained in the same way as described in Section 2. 

One may ask whether it is useful to implement a method to compute grids on minimal surfaces in a 
3D multi-block grid generator code. The answer is yes. Minimal surfaces may be used to define the 
geometry and grid for a block-face of which only the four face-edges are given. It is also possible 
to apply minimal surface grid generation when a grid must be generated in a block-face with four 
face-edges lying in a plane. Then the minimal surface is a plane surface bounded by the four edges. 
The grids in the 2D domains depicted in Figs. 5, 6, 7 were generated in this way and are in fact grids 
on minimal surfaces. 

An example of a grid on a characteristic minimal surface is shown in Fig. 8. This is a so-called 
square Scherck surface [8]. Fig. 9 illustrates what happens when the prescribed boundary grid 
point distribution is changed. This figure clearly shows that the shape of the minimal surface is 
independent of the prescribed boundary grid point distribution. 

4 Surface Grid Generation on Parametrized Surfaces 

In this section we develop a method to generate a grid on a parametrized surface which is indepen- 
dent of the parametrization. A generated grid only depends on the shape of the surface and the 
prescribed boundary grid point distribution at the four edges of the surface. 

Consider a bounded surface S with a prescribed geometrical shape in three dimensional physical 
space with Cartesian coordinates z = (x,y,z) T . Assume that S is parametrized by a differentiable 
one-to-one mapping 

x:V uv »S, (42) 

where V uv is the unit square in two dimensional space with Cartesian coordinates u = (u, v) T . 
Define the four edges E\, E 2 , E 3 , E\ of surface S by 

• u = 0 at edge E\ and u = 1 at edge E 2 , 
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Figure 3: Transformation from computational (£,77) space to a parametrized surface S in Cartesian 
(#, 3/, z) space. 

• v = 0 at edge £3 and v = 1 at edge £4. 

Thus (£i,£ 2 ) and (£3, £4) are the two pairs of opposite edges of surface S as shown in Fig. 3. 
Introduce the parameter space V s t as the unit square in a two dimensional space with Cartesian 
coordinates s = (s,^. Again require that the parameters s and t obey the boundary equations 
specified in Eqs.(l),(2),(3),(4). Furthermore, require that As = 0 and At — 0 where A is the 
Laplace- Beltrami operator for surfaces. Hence the parameters s and t obey 

T J o} 2 Sy^ -(“ o} 2 s u T Ja 22 Sy^ — 0, (43) 

(ja u t n + Ja l2 t v S j + ( Ja l2 t u + Ja 22 t v ^ = 0, (44) 

where a %3 are the contravariant tensor components and J 2 is defined as the determinant of the 
covariant metric tensor. The contravariant tensor components a 13 are related to the covariant 
tensor components according to Eq.(15). The covariant metric tensor components are defined 
by Eq.(12), where the two covariant base vectors are now given by 

a\ = x u , a 2 = x v . (45) 

The coefficients Ja 11 , Ja 12 and Ja 22 in Eqs.(43),(44) are thus functions of u and u, and Eqs.(43),(44) 
are therefore two uncoupled second-order linear partial differential equations for the functions 
s — s(u, n) and t = /(u, u). 

Each boundary point of surface S has a unique (s, t) parameter value at dV s t and a unique (u, n) pa- 
rameter value at dV uv . Thus each (u, u) parameter value at dV uv has also a unique (5, t) parameter 
value at dV s t • Thus the functions s and t are prescribed at the boundary of V U v Hence, Eq.(43) 
together with the Dirichlet boundary conditions for s can be used to compute s = s(u,v) and 
Eq.(44) together with the Dirichlet boundary conditions for t can be used to compute t = t(u,v). 
Only two linear partial differential equations have to be solved to define these mappings. These 
two mappings are compactly written as s : V uv »— ► V st . Note that s : V uv V s t is a differentiable 
one-to-one mapping so that the inverse mapping u : V s t »— ► V uv also exists. 

Thus the composite mapping x : V s t ^ S , defined as x = f(u(^)) also exists and is differentiable 
and one-to-one. Note that this mapping x : V s t ► S only depends on the shape of surface S and 
is independent of the original parametrization x : V uv S. The mapping x : V st ^ S may thus 
be considered as a property of surface S and defines a new unique parametrization of S. 
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Consider a prescribed boundary grid point distribution at the four edges E\, E 2 , E 3, £* 4 . The bound- 
ary grid point distribution can be defined as a mapping x : dC dS where C is the computational 
space defined as the unit square in a two dimensional space with Cartesian coordinates £ = (f , j]) T . 
Because x : dC »— ► dS is prescribed and x : dV s t »-► $5 is defined as described above, it follows that 
s : dC »— ► <9P S f is also defined. 

In exactly the same way as for the two dimensional case, the mapping ?:Ch V st is now defined by 
the algebraic straight line transformation defined by Eqs.(9),( 10). The composition of the mapping 
s : C 1 ► and the mapping £ : V s t 5 defines x : C 5 and describes the interior grid 
point distribution on surface S. Note that this composite mapping will also be differentiable and 
one-to-one. 

Grid orthogonality at boundaries can be obtained by the same procedure as described in Section 2 
for 2D domains. 

Fig. 10 shows an irregularly distributed control point mesh on a smooth surface. The surface is de- 
fined as z — ^tanh( 15(^-(x- l) 2 — (y- l) 2 )), (x, y) E [0, l] 2 . In section 5.2 is described how bicubic 
Hermite interpolation is used to define the mapping x : V uv ► S . The parametrization depends on 
the control point distribution. Fig. 11 shows an elliptic grid. Equidistributed boundary grid points 
are used as Dirichlet boundary condition. This figure clearly demonstrates that the interior surface 
grid only depends on the shape of the surface and is independent of the parametrization. 

Less academic surface grids are shown in Fig. 13. 

5 Surface modeling 

Consider a discrete surface defined as a two-dimensional array of control points. An interpolated 
surface is obtained by bicubic Hermite interpolation. Non-linear averaging formulas for the tangent 
and twist vectors are used to prevent spurious oscillations. The interpolation method is explained 
first for curves and then extended to surfaces. 


5.1 Piecewise cubic Hermite interpolation for curves 


Consider a set of control points {x, = (x,y,z)f | i = 0. . .N}. We wish to construct a smooth 
C l curve x:uE[0, 1]h 1Z 3 which is passing through the set of control points with a geometrical 
shape as one would intuitively expect. Furthermore, spurious oscillations should be prevented. For 
this reason, cubic spline interpolation is not safe. Instead, piecewise cubic Hermite interpolation is 
applied. The extra freedom to model the tangent vectors is used to prevent unwanted oscillations. 
The parameter u is defined as normalized arc length. 

Compute the distance between succeeding control points: 


d{ =|| Xi - x t _ ! || , i - 1 . . . N. 


(46) 


Define the length of the curve by 


N 


£ = £4 

t=l 


(47) 


627 



and the normalized distances as 


( 48 ) 


di =■ di/L , i = 1 . . .N. 

Define the knot sequence {u, | i = 0 . . . AT) by uq = 0 and 

u, = u t _i + di , i = 1 . . . N. (49) 

Hence, 0 = uo < U\ < . . . < u n = 1. Patch i is defined as the interval [«,•_!,«,]. In patch i, we 
relate a local variable a G [0, 1] to the global variable u by 

u = «,_! + a(ui - «,•_!) = ««•_! + adi. (50) 

For the moment, suppose that the tangent vectors {x Ui = ^(u t ),i = 0...N } are known. How 
these tangent vectors are modeled is shown below. 


The shape of the curve at patch i is then defined as 

£(a) = fj_!^ 0 (a) + XiHi(a) + d i x Ut _ 1 H 2 {a) + d,x Ui // 3 (a), 
where Hq , H i, H 2 , H 3 are cubic Hermite interpolation functions defined as 

H 0 (a) = (1 + 2a)( 1 — a) 2 , 

H\(a) = (3 — 2 a)a 2 , 

H 2 (a) = o(l -«) 2 , 

H 3 (a) = (a-l)a 2 , 

with 0 < a < 1. 


(51) 


(52) 


It can be easily verified that £(«.-> = jf( = x Ui , so that the piecewise cubic curve x(u) is 
indeed C 1 . 


The tangent vectors {x Ut , i = 0. . .N} are computed as follows. Define the chord vectors 

_ X{ — £V_i 

= j , * = 1 - . - N. ( 53 ) 

Note that || ||= L. The tangent vectors at the interior knots i = 1 . . .N — 1 are modeled as 

i _ i T u ^ ( 1 C|) j i — l . . .N 1, ( 54 ) 

with 


c* 


— 1 


df 


, i = I...N - 1. 


(55) 


ii ~ *i-i II 2 + || *.+i - *« || 2 df + df +1 

If || x,- — Xj_ j ||<|| x^+j — x,- || then ss 0 and x Ui « This implies that high curvature 

will be restricted to small patches which is a behaviour as one would intuitively expect. Spurious 
oscillations are also prevented. 


Quadratic end conditions are used to compute the end tangent vectors x UQ and x UN . The quadratic 
end conditions require that the cubic polynomial function x{a) is a quadratic function of a in patch 
1 and in patch N . It is easily verified that this implies that 


= 2x 




‘'ll! 


%u N ~ 


(56) 


Fig. 12 illustrates that cubic Hermite interpolation prevents spurious oscillations, in contrast to 
cubic spline interpolation. 
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5.2 Piecewise bicubic Hermite interpolation for surfaces 

Consider a set of control points = (x, y , z)Jj \ i — 0 . . . N, j — 0 . . . M }. We wish to construct 
a smooth C A surface x : (u,u) £ [0, l] 2 »-► 1Z 3 which is passing through the set of control points 
with a geometrical shape as one would intuitively expect. As for curves, spurious oscillations should 
be prevented. For this reason, bicubic spline interpolation is not safe. Instead, piecewise bicubic 
Hermite interpolation is applied. The extra freedom to model the tangent vectors and twist vectors 
is used to prevent unwanted oscillations. 

Consider a row of points | i = 0 ...N} with j £ {0...M} fixed. This row of points is 

considered as a discrete curve and it is therefore possible to compute a knot sequence | i = 

0...N} in exactly the same way as described in the previous section. In the same way, consider 
a column of points {x hJ | j = 0...M} with i £ {0...7V} fixed, and compute the knot sequence 

I j = 

To construct a smooth surface, one knot sequence is needed for all rows and all columns. These 
two knot sequences are obtained by averaging: 

1 M j N 

Ui = a7TT^ U,j ’ * = 0...N, V, = » 3 = (57) 

m -t i i t a . =0 

Patch (i, j) is defined as the rectangle u x ] x Vj]. In patch (i, j) we relate local variables 

(a,/3) £ [0, l] 2 to the global variables (u,u) by 

u = Ui-i + ad^ , v - 1 + fid", (58) 

with d™ = U{ — Ui - 1 and d v - ~ vj — Vj-\. 

For the moment, suppose that the tangent vectors x Ut} = | ^(ui,Vj) , x Vt = ^(u ty Vj) , and twist 
vectors x UVtj = ( U{ , Vj) are known for all knots (i, j). How these tangent and twist vectors are 

modeled is shown below. 

The shape of the surface at patch (i, j) is then defined as 


x(a,0) = ( Ho(a), H,(a), H 2 (a), H 3 (a))M% 


H 0 (P) \ 

am 

urn 

am 


where the matrix M,^ is defined by 



1 

7 

n? 


d-x v , . 



1 


iVV J 

d j £ v tJ 


T 

7 

3-~ 

^3 

//V r 
u » 

d u d v x 

d ld)x uv> _^ 


fjV'X 

u t 

d u x 

a % X Ui,J 

d u d v x 

d t d )x UVl] / 


(59) 


(60) 


From these definitions, it can be easily verified that the piecewise bicubic surface f(u,t;) is C 1 . 
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The tangent vectors x Ut are computed as follows ( the tangent vectors x Vt j are computed in the 
same way ). Define the chord vectors 


-* x i,j x i—\ jj 

“ d? 


, i = 1 . . .N , j = 0 . . .M, 


and use the same non-linear averaging formula as used for curves, thus 


(61) 


X u x%3 — "I ~ *^ u *+ £ j ( ^ ^ — h . . iV 1?J — 0... M. 


with 


Ci.j 


x tJ 


“I" || X i+\J x iJ 


, i = 1 . . . N — \ , j — 0 . . . A/. 


(62) 

(63) 


Quadratic end conditions are used again to compute the end tangent vectors. 


A modification of Adini’s method [9] is used to compute the twist vectors. Consider patch (i, j) with 
local variables (a,/?). Assume that the tangent vectors x u , x v are known at the four corner points of 
the patch. Introduce the abbreviate notation x 0 o — Xjo = x i,j- l, x o\ = x i- ij, Xu = 

Use Eqs.(59),(60) to find £ a (0, 0) = d“x Ut _ x , f a (l , 0) = c^x*. , x a (0, 1 ) = d“x u ,_, ; , x a ( 1 , 1 ) = 
d?x Utj , ^(0,0) = djx„,_, x>(l,0) = d V jX Vil _ % , f^(0, 1) = d)x Vi _^ } , xp(\, 1) = Compute 

the boundary curves of the patch by cubic Hermite interpolation. Thus, for example 


£(a,0) = foo^o(a) + x X0 H\(a) + x a (0,0)tf 2 (a) + x a (l, 0)tf 3 (a), (64) 


and the other three boundary curves are computed by similar formulas. Define the shape of the 
surface patch as a bilinearly blended Coon’s patch 


x(a,P) = (1 -a, a) 


x(0,/3) 

x(l,/J) 


+ (1 ~P,P) 


x(a, 0) 
x(a, 1) 


(1 




x OO x O\ 
X\o Xu 



(65) 


Compute the corner twist vectors x a p(0, 0), z a ^(l, 0), x a /?(0, 1), x a p{\, 1) from Eq.(65). Use Eq.(58) 
to find the corresponding twist vectors w.r.t. the global variables (u,u): x uv% _ V j _^ — x a ^(0,0 )/d^dj 
etc.. 


Thus at the four corners of each patch, an estimation is found for the twist vector x uv . Consider 
an interior knot (z,j). Then four estimations for the twist vector x UVt are found in respectively 
patches +1) and (i -f 1,J + 1)- Write those estimations as respectively 

^uvi % uv %uv E - A similar averaging procedure as applied for tangent vectors is used to define a 
unique value x UVxj : 


fSW a2 , ?SE a2 I ~NW a 2 , t NE A? 

x uv * x uv ^i+l j * ' x uv ^t+lj+1 

A h + A f+hj + A h + 1 + ^i+ij+i 


where Atj is the area of patch (i,j) defined as 


— 0.5 || ( x ij A 


( 66 ) 


(67) 


This non-linear averaging procedure guarantees that large changes in twist will be restricted to 
small patches. At a boundary knot, there are only two estimations for the twist vector. It is 
evident how the averaging procedure must be applied in that case. At the four corner knots, only 
one estimation is available and averaging is thus not needed. 


Figs. 10, 11 illustrate that bicubic Hermite interpolation gives a smooth surface shape, even if the 
control points are irregularly distributed. 
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Concluding remarks 


An elliptic grid generation method is presented to generate boundary conforming grids in domains 
in 2D physical space and on minimal surfaces and parametrized surfaces in 3D physical space. The 
method is based on the use of composite mappings and produce excellent grids in the sense of 
smoothness, grid point distribution and regularity. 

The described elliptic grid generation method has been implemented into NLR’s multi-block grid 
generation code ENGRID and is extensively used for the generation of boundary conforming Navier- 
Stokes grids in block-faces with complex shapes. 

Concerning surface modeling, it is shown that bicubic Hermite interpolation is an excellent method 
to define interpolated surfaces. 
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Figure 4: The NLR ENFLOW system 



Figure 5: Elliptic grid with grid orthogonality at the lower boundary of the domain. 
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Figure 6: Elliptic grid with boundary layer and orthogonality. 



Figure 7: Detail of elliptic grid at convex part of the boundary. 
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Figure 8: Minimal surface grid. Surface is a square Scherck surface. 



Figure 9: Minimal surface grid. Shape of surface is independent of the boundary grid point distri- 
bution. 
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Figure 10: Irregularly distributed control point mesh on a smooth surface. 



Figure 11: Elliptic grid on the surface. Grid is independent of the parametrization. 
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Figure 12: Comparison of cubic Hermite and cubic spline interpolation. 



Figure 13: Surface grids of a wing-body-pylon-nacelle configuration. 
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